library("poppr")
library("mmod")
library("magrittr")
library("pegas")
library("hierfstat")
#Cargar el archivo guardado en formato delimitado por ","
data.maiz <- read.csv2("maiz.csv")
#Crear el objeto loci
alelos.maiz <- alleles2loci(data.maiz, ploidy = 2, rownames = 1, population = 2)
#Crear el objeto genind
magen <- loci2genind(alelos.maiz, ploidy = 2, na.alleles = "NA")
#Asignar estratos desde información de otra tabla datos
estratos <- read.csv("MaGene_Strata.csv", sep =";", header = TRUE) #importante revisar que la tabla de datos tenga la información ordenada según los individuos
strata(magen) <- estratos #estamos asignando la tabla de datos que importamos directamente
nameStrata(magen) <- ~region/sitio/color #para que lo haga de manera anidada Laboratorio 2 - Genética de Poblaciones
Estimativas de diversidad genética
Objetivo: aprender a manejar datos genéticos con estratificación y calcular las estimativas básicas de diversidad genética.
Para los ejemplos, se va a utilizar la base de datos de maíz criollo que consiste en 224 muestras de provenientes de las regiones Brunca y Chorotega de Costa Rica, analizadas con 20 microsatélites (SSR), 10 son dinucleotídicos (bnlg) y 10 son tetranucleotídicos (phi). Dentro de cada región, se establecieron sitios de muestreo (Fig.3) y dentro de cada sitio se muestrearon accesiones las cuales consisten en cultivares o milpas de un tipo de maíz criollo de diferente color: morado, congo (morado oscuro), blanco o amarillo. El color se considera un factor de estructuración, por lo que puede analizarse como un estrato más. Los individuos corresponden a mazorcas dentro de cada color, accesión o sitio. Los análisis a continuación deberán realizarse en alguno de los estratos de interés: región, sitio o color de mazorca.
Fig 3. Mapa de Costa Rica con las regiones Chorotega (al norte) y Brunca (al sur) resaltadas en gris. En las regiones se muestran de color los sitios de colecta. Dentro de la Región Chorotega cada círculo corresponde a una accesión en cada sitio de colecta; en la Región Brunca, las accesiones se representan con cuadraros
1. Importar datos y seleccionar poblaciones
Se procede a importar los datos, a crear el objeto genind y a asignar la estratificación como se cubrió en la práctica 1. Recuerden que las librerías para análisis genéticos deben cargarse con el comando library(). Refiérase a la práctica 1 para la lista de librerías. Recuerde que debe definir el directorio de trabajo, con el comando setwd() o utilizando explorador de archivos de Rstudio.
Cuando se desea hacer los análisis para un estrato en específico (sitio o color) se debe indicar en R en cuál de estos niveles se desea trabajar. El comando setPop(x) <- value que vimos en el laboratorio anterior, permite cambiar el estrato que se quiere utilizar. En value, se indica una fórmula, que consiste en una virgulilla “~” seguida del nombre del estrato donde R calculará las estadísticas de diversidad. Por ejemplo, si se quisiera analizar la diversidad por región socioeconómica value sería ~region. Si se quisiera analizar por región y sitio de colecta se indica ~region/sitio. El objeto magen tiene las columnas con 3 estratos, por lo que para cada análisis debe hacerse esta asignación de estrato. A continuación usamos el comando para definir la región (i.e., Brunca o Chorotega) como el estrato de interés. Es decir, las estadísticas de diversidad se calcularán por separado para la región Brunca y para la región Chorotega. En este análisis, el sitio y el color NO se toman en cuenta, porque no los estamos incluyendo en la fórmula:
setPop(magen) <- ~region 2. Información por locus
2.1. Prueba del Equilibrio Hardy-Weinberg
Usando el comando hw.test de la librería pegas, podemos probar la hipótesis de que las frecuencias genotípicas de una serie de loci siguen el equilibrio de Hardy-Weinberg. Lo que realiza este comando es una prueba de chi-cuadrado entre las frecuencias genotípicas observadas y las esperadas para cada SSR:
hw.test(alelos.maiz) # Incluye todos los genotipos de un SSR chi^2 df Pr(chi^2 >) Pr.exact
bnlg1046 2111.4313 231 0 0
bnlg1191 1118.6373 300 0 0
bnlg1194 2345.4324 703 0 0
bnlg1265 1140.3796 406 0 0
bnlg1449 1199.2219 136 0 0
bnlg1523 2250.3223 595 0 0
bnlg1740 2131.9848 378 0 0
bnlg1890 1775.9587 528 0 0
bnlg1917 2419.4205 351 0 0
bnlg2305 532.7068 136 0 0
phi015 448.4887 55 0 0
phi031 254.2561 10 0 0
phi064 758.2628 136 0 0
phi072 440.5250 55 0 0
phi078 930.8504 55 0 0
phi089 481.4608 10 0 0
phi093 1064.8512 55 0 0
phi109188 340.9622 45 0 0
phi116 419.5332 36 0 0
phi96100 465.9575 45 0 0hw.test(magen) # No incluye homocigotos, observar que los grados de libertad (df) son menores que cuando se prueba en alelos.maiz chi^2 df Pr(chi^2 >) Pr.exact
bnlg1046 1592.52013 210 0.000000e+00 0.000
bnlg1191 826.74072 276 0.000000e+00 0.000
bnlg1194 1865.72402 666 0.000000e+00 0.000
bnlg1265 850.92391 378 0.000000e+00 0.000
bnlg1449 814.13613 120 0.000000e+00 0.000
bnlg1523 1827.30854 561 0.000000e+00 0.000
bnlg1740 1635.41552 351 0.000000e+00 0.000
bnlg1890 1461.88968 496 0.000000e+00 0.000
bnlg1917 1803.38114 325 0.000000e+00 0.000
bnlg2305 272.87474 120 6.405987e-14 0.000
phi015 213.46473 45 0.000000e+00 0.000
phi031 27.68975 6 1.074820e-04 0.002
phi064 448.39915 120 0.000000e+00 0.000
phi072 203.95877 45 0.000000e+00 0.000
phi078 624.80525 45 0.000000e+00 0.000
phi089 234.47321 6 0.000000e+00 0.000
phi093 758.26764 45 0.000000e+00 0.000
phi109188 111.74069 36 1.096105e-09 0.000
phi116 177.20194 28 0.000000e+00 0.000
phi96100 223.59467 36 0.000000e+00 0.000La columna Pr(chi^2 >) corresponde a los valores de p-value para cada SSR. Es decir que Pr(chi^2 >) sería la probabilidad de obtener un valor de chi-cuadrado mayor que el p-value bajo la hipótesis nula (Ho=La población se encuentra en equilibrio Hardy-Weinberg). Si Pr(>Chisq) es mayor que el p-value, se rechaza Ho, es decir, se concluiría que la población no está en HWE. Para los resultados de esta tabla, vemos que Pr(>Chisq) = 0, por lo que se acepta Ho, asumiendo que la población está en HWE.
2.2. Estimativas con poppr
Con el comando locus_table(x, population = "ALL")de poppr se obtienen las estimativas de diversidad para cada marcador, ya sea para todas las poblaciones (por default o indicando population = "ALL") y para alguna población en específico (population = "Brunca"). El resultado es una tabla de cuatro columnas donde se indica el número de alelos observados (allele), índice de diversidad (default es el Simpson 1-D) y la heterocigosidad esperada (Hexp). El cálculo para la heterocigosidad esperada es \frac{n}{n-1} \sum{1-p_i} donde p es la frecuencia de alelos en un locus dado y n es el número de alelos observados (Nei, 1978). También se obtiene el “Eveness” que es el índice de Simspon expresado como un porcentaje de la diversidad máxima. Es decir, entre más cercano a 1, más diversidad hay.
locus_table(magen) #por default lo hace para todas las poblaciones en conjunto
allele = Number of observed alleles
1-D = Simpson index
Hexp = Nei's 1978 gene diversity
------------------------------------------ summary
locus allele 1-D Hexp Evenness
bnlg1046 21.00 0.87 0.87 0.61
bnlg1191 24.00 0.94 0.94 0.85
bnlg1194 37.00 0.96 0.96 0.84
bnlg1265 28.00 0.94 0.94 0.81
bnlg1449 16.00 0.89 0.89 0.82
bnlg1523 34.00 0.90 0.90 0.56
bnlg1740 27.00 0.90 0.90 0.67
bnlg1890 32.00 0.87 0.88 0.49
bnlg1917 26.00 0.81 0.81 0.46
bnlg2305 16.00 0.88 0.88 0.84
phi015 10.00 0.76 0.76 0.68
phi031 4.00 0.66 0.66 0.84
phi064 16.00 0.90 0.90 0.83
phi072 10.00 0.73 0.74 0.64
phi078 10.00 0.75 0.75 0.75
phi089 4.00 0.49 0.49 0.79
phi093 10.00 0.79 0.79 0.70
phi109188 9.00 0.78 0.78 0.85
phi116 8.00 0.80 0.80 0.88
phi96100 9.00 0.77 0.77 0.76
mean 17.55 0.82 0.82 0.73En el cuadro anterior se observa la diversidad por locus. Vemos que el número de alelos varía entre 4 y 37, con un promedio de 17.55 alelos. Estos valores son bastante altos, sin embargo, para microsatélites se esperan muchos alelos. En relación a la heterocigosidad esperada (Hexp), se observa que en promedio He = 0.82. Es decir, el 82% de los individuos son heterocigotos. Para microsatélites (SSR) se esperan valores de heterocigosidad superiores a 0.7, ya que este marcador tiene un leve sesgo de búsqueda, en el cual se seleccionan microsatélites que son altamente polimórficos. Los índices de Simpson y Evenness, rara vez se usan en publicaciones y se presentan aquí por factores históricos.
3. Estimativas de diversidad por población
Para obtener las estimativas de diversidad genética, existen varias librerías y funciones que se pueden utilizar. Entre estas poppr, adegenet, hierfstat, pegas. En esta práctica veremos algunas de las funciones dentro de estas librerías.
3.1. Función summary
El comando summary(x) se utiliza en muchos otros análisis. En nuestro caso va a mostrar: el número de individuos, número y tamaño de cada población, el número de alelos por locus, el número de alelos por población, el total de valores perdidos y las heterocigosidades esperadas y observadas por locus. Se recalca que si de desean estas estimativas para un estrato en específico se utiliza el comando setPop(x) <- value y luego el comando summary(x).
setPop(magen) <- ~region
summary(magen)
// Number of individuals: 224
// Group sizes: 81 143
// Number of alleles per locus: 21 24 37 28 16 34 27 32 26 16 10 4 16 10 10 4 10 9 8 9
// Number of alleles per group: 274 323
// Percentage of missing data: 10.27 %
// Observed heterozygosity: 0.34 0.66 0.57 0.71 0.42 0.59 0.47 0.55 0.34 0.61 0.66 0.59 0.58 0.57 0.31 0.33 0.4 0.61 0.51 0.65
// Expected heterozygosity: 0.87 0.94 0.96 0.94 0.89 0.9 0.9 0.87 0.81 0.88 0.76 0.66 0.9 0.73 0.75 0.49 0.79 0.78 0.8 0.77Uno puede calcular el promedio de la heterocigosidad observa y esperada para todas las poblaciones juntas, si calcula un promedio. Sin embargo, veremos más adelante, que esto puede conllevar problemas por un efecto Wahlund.
mean( summary(magen)$Hexp )[1] 0.8186826sd( summary(magen)$Hexp )[1] 0.11144473.2. Función allelic.richness
Función que pertenece a la librería hierfstat. Estima la riqueza alélica, los recuentos alélicos rarificados, por locus y población. Esta estadística es importante cuando se trabaja con poblaciones que difieren mucho en el tamaño de muestra. En estos casos, estadísticas como el número de alelos pueden estar fuertemente sesgadas en muestras con poco número de individuos. Por ello se usa la rarefacción. Estadísticas como la heterocigosidad, son un poco más robustas a la diferencia en tamaños de muestra. Para calcular el número de alelos mediante rarefacción, se utiliza de la forma allelic.richness(data, diploid = TRUE). Los argumentos consisten en un objeto genind (data) y si son datos diploides (default) o no. Como resultado devuelve min.all que es el mínimo tamaño poblacional utilizado para la rarefacción y una tabla con tantas filas como loci y tantas columnas como poblaciones que contienen la cuenta de alelos rarificados.
riqueza <- allelic.richness (magen)
class(riqueza)# corroborar que es un objeto tipo lista [1] "list"riqueza $min.all
[1] 134
$Ar
BRUNCA CHOROTEGA
bnlg1046 13.943533 19.048891
bnlg1191 22.316156 20.631796
bnlg1194 30.695082 30.523445
bnlg1265 23.370204 23.346954
bnlg1449 11.910705 14.381959
bnlg1523 23.396699 26.604865
bnlg1740 18.795197 21.941836
bnlg1890 22.264797 23.339412
bnlg1917 16.000000 20.047552
bnlg2305 13.879566 11.950832
phi015 8.985205 9.825632
phi031 4.000000 3.997477
phi064 11.930531 15.131251
phi072 7.856185 9.612051
phi078 8.775838 9.166508
phi089 3.948899 2.999853
phi093 8.000000 9.515764
phi109188 6.825423 6.977148
phi116 5.999999 7.297079
phi96100 6.948899 8.493664Aquí se presenta el número de alelos rarificado por locus, en dos columnas que representan Brunca y Chorotega (siempre R acomoda los factores en orden alfabético). Para obtener una estimativa por región, se usa el comando ColMeans(riqueza$Ar). Con este comando se observa que la riqueza alélica para Brunca es Ar = 13.492 y para Chorotega es Ar = 14.742. Hay mayor riqueza alélica en Chorotega, sin embargo, las estadísticas son bastante similares ya que entre ellas no hay una gran diferencia en el tamaño de población.
También se puede utilizar la funcion apply para estimar tanto el promedio como la desviación estandard de la riqueza alélica. Esta función requiere de tres argumentos. El primero es un marco de datos, que en este caso es la riqueza alélica por locus y región que se encuentra en riqueza$Ar, luego se debe indicar si vamos a aplicar una función a las filas (1) o las columnas (2). Y finalmente se incluye la función a utilizar, en este caso son dos funciones: mean y sd.
apply(riqueza$Ar, 2, mean) BRUNCA CHOROTEGA
13.49215 14.74170 apply(riqueza$Ar, 2, sd) BRUNCA CHOROTEGA
7.661231 7.933973 3.3. Función basic.stats
Permite obtener estimativas básicas de diversidad genética (Cuadro 1). Se utiliza: basic.stats(data, diploid=TRUE, digits=4). Se le debe indicar con cuál objeto genind se trabajará en data, si los datos son diploides (default) y el número de dígitos de las estimativas. Los resultados pueden verse como una lista de las estimativas por locus (filas) y poblaciones (columnas) o como un cuadro de datos (data frame) de las estimativas (columnas) por cada locus (filas). Como se ha mencionado, para acceder a algún elemento de la lista se usa el nombre del objeto donde se guardaron las estadísticas básicas, el símbolo $ y luego la abreviación de la estimativa (Cuadro 1).
Cuadro 1. Descripción de estimativas calculadas por la función basic.stats. En negrita se resaltan los de interés.
| Abreviación | Estimativa |
|---|---|
| Ho | Heterocigosidad observada |
| Hs | Heterocigosidad esperada dentro de la población H_S = 1 - \sum{p_i} donde p_i es la frecuencia del i-ésimo alelo. Esta medida es equivalente a la frecuencia de los heterocigotos bajo las expectativas de Hardy-Weinberg |
| Ht | Diversidad genética total. H_T = 1 - \sum{p_i} , donde p_i es la frecuencia media del i-ésimo alelo en la muestra agrupada obviando los estratos. |
| Dst | Diversidad genética entre las poblaciones |
| Htp | Diversidad genética general corregida por el tamaño de muestra |
| Dstp | Diversidad genética entre las muestras corregida por el tamaño de muestra |
| Fst | Índice de Fijación: mide el déficit de heterocigotas de las subpoblaciones en relación con la población total. Valores entre 0 y 1 |
| Fstp | F_{ST} corregido por el tamaño de muestra |
| Fis | Endogamia, mide la reducción en la heterocigosidad de los individuos dentro de las sub-poblaciones valores entre -1 y 1. |
| Dest | Medida diferenciación genética según Jost (2008) |
Con los siguientes comandos se obtienen las estadísticas básicas generales:
bs <- basic.stats(magen, digits= 4)
bs $perloc
Ho Hs Ht Dst Htp Dstp Fst Fstp Fis Dest
bnlg1046 0.3287 0.8720 0.8748 0.0028 0.8776 0.0056 0.0032 0.0064 0.6230 0.0437
bnlg1191 0.6587 0.9332 0.9381 0.0049 0.9430 0.0098 0.0052 0.0104 0.2941 0.1470
bnlg1194 0.5796 0.9578 0.9623 0.0044 0.9667 0.0089 0.0046 0.0092 0.3949 0.2103
bnlg1265 0.7051 0.9394 0.9447 0.0053 0.9499 0.0105 0.0056 0.0111 0.2494 0.1740
bnlg1449 0.4111 0.8739 0.8906 0.0167 0.9074 0.0334 0.0188 0.0368 0.5296 0.2650
bnlg1523 0.5942 0.8948 0.9016 0.0069 0.9085 0.0137 0.0076 0.0151 0.3360 0.1304
bnlg1740 0.4590 0.8796 0.8957 0.0161 0.9119 0.0322 0.0180 0.0354 0.4782 0.2678
bnlg1890 0.5512 0.8689 0.8768 0.0079 0.8846 0.0157 0.0090 0.0178 0.3657 0.1198
bnlg1917 0.3266 0.8034 0.8099 0.0066 0.8165 0.0132 0.0081 0.0161 0.5934 0.0670
bnlg2305 0.5942 0.8798 0.8861 0.0062 0.8923 0.0124 0.0070 0.0139 0.3247 0.1035
phi015 0.6574 0.7456 0.7643 0.0187 0.7830 0.0374 0.0245 0.0478 0.1183 0.1470
phi031 0.5895 0.6682 0.6703 0.0021 0.6724 0.0042 0.0031 0.0062 0.1179 0.0126
phi064 0.5620 0.8935 0.9004 0.0069 0.9073 0.0138 0.0077 0.0152 0.3710 0.1294
phi072 0.5636 0.7287 0.7315 0.0029 0.7344 0.0057 0.0039 0.0078 0.2266 0.0210
phi078 0.3156 0.7263 0.7460 0.0197 0.7657 0.0394 0.0264 0.0514 0.5654 0.1439
phi089 0.3208 0.4898 0.4955 0.0056 0.5011 0.0113 0.0114 0.0225 0.3452 0.0221
phi093 0.4180 0.7786 0.7952 0.0166 0.8118 0.0333 0.0209 0.0410 0.4631 0.1502
phi109188 0.5902 0.7572 0.7688 0.0116 0.7804 0.0233 0.0151 0.0298 0.2205 0.0958
phi116 0.5001 0.7782 0.7911 0.0129 0.8040 0.0258 0.0163 0.0320 0.3574 0.1161
phi96100 0.6617 0.7542 0.7619 0.0077 0.7696 0.0154 0.0101 0.0200 0.1226 0.0625
$overall
Ho Hs Ht Dst Htp Dstp Fst Fstp Fis Dest
0.5194 0.8112 0.8203 0.0091 0.8294 0.0182 0.0111 0.0220 0.3597 0.0966 bs$Ho #obtener la Heterocigosidad observada BRUNCA CHOROTEGA
bnlg1046 0.2676 0.3898
bnlg1191 0.6456 0.6719
bnlg1194 0.6027 0.5565
bnlg1265 0.6795 0.7308
bnlg1449 0.3836 0.4386
bnlg1523 0.5946 0.5938
bnlg1740 0.4058 0.5122
bnlg1890 0.5385 0.5639
bnlg1917 0.2687 0.3846
bnlg2305 0.5270 0.6613
phi015 0.6579 0.6569
phi031 0.6053 0.5736
phi064 0.4861 0.6379
phi072 0.5256 0.6015
phi078 0.3467 0.2846
phi089 0.2625 0.3790
phi093 0.4868 0.3492
phi109188 0.5063 0.6741
phi116 0.4533 0.5469
phi96100 0.7250 0.5984bs$Fis #obtener solamente el índice de endogamia BRUNCA CHOROTEGA
bnlg1046 0.6996 0.5432
bnlg1191 0.3040 0.2844
bnlg1194 0.3701 0.4196
bnlg1265 0.2751 0.2238
bnlg1449 0.5612 0.4981
bnlg1523 0.3416 0.3302
bnlg1740 0.5195 0.4400
bnlg1890 0.3771 0.3544
bnlg1917 0.6577 0.5320
bnlg2305 0.4049 0.2436
phi015 0.1116 0.1248
phi031 0.1258 0.1095
phi064 0.4603 0.2804
phi072 0.2592 0.1956
phi078 0.5039 0.6224
phi089 0.4913 0.1830
phi093 0.3865 0.5426
phi109188 0.2940 0.1543
phi116 0.3897 0.3278
phi96100 -0.0188 0.2486bs$overall #obtener el promedio de cada estimativa para todas las poblaciones Ho Hs Ht Dst Htp Dstp Fst Fstp Fis Dest
0.5194 0.8112 0.8203 0.0091 0.8294 0.0182 0.0111 0.0220 0.3597 0.0966 write.csv(bs$perloc,"basic.stats.csv")#exportar los dato con todas las estimativasAl igual que en resultados anteriores podemos interpretar varios resultados. Esta tabla nos muestra que para la muestra total, hay una heterocigosidad observada de H_O = 0.519. Es decir, cerca del 60% de la población es heterocigota. La probabilidad de encontrar un heterocigoto según HWE es de 83% (Hs=0.8112) por lo que tenemos un déficit de heterocigotos en la población total. Esto veremos más adelante se puede deber a un efecto Wahlund porque estamos mezclando dos regiones muy diferentes, la Brunca y la Chorotega. Lo más correcto es tener estimativas por región. Se pueden utilizar las herramientas de selección de listas, matrices y de tablas de datos de R para obtener información específica por población:
fis <- bs$Fis # extrae solo el Fis y lo guarda en un objeto
Brunca_fis<- bs$Fis[ ,1] # utiliza la selección de columnas para tomar solo los datos de brunca
Brunca_fis bnlg1046 bnlg1191 bnlg1194 bnlg1265 bnlg1449 bnlg1523 bnlg1740 bnlg1890
0.6996 0.3040 0.3701 0.2751 0.5612 0.3416 0.5195 0.3771
bnlg1917 bnlg2305 phi015 phi031 phi064 phi072 phi078 phi089
0.6577 0.4049 0.1116 0.1258 0.4603 0.2592 0.5039 0.4913
phi093 phi109188 phi116 phi96100
0.3865 0.2940 0.3897 -0.0188 class(Brunca_fis) [1] "numeric"bfis <- as.data.frame(Brunca_fis) #para guardarlo como tabla. Normalmente, como interesa analizar las estimativas por población y no tanto por marcador, se trabaja con promedios de las distintas estimativas. Para esto se utiliza el comando apply(x, MARGIN, FUN), que aplica una función (FUN) a las filas o columnas (MARGIN) de una matriz o lista de valores. Por ejemplo, MARGIN =1 indica que la función se aplicará a filas, mientras que MARGIN=2 indica que la función FUN se aplicará a las columnas. El argumento FUN permite especificar la función que se desea calcular, por ejemplo, promedio (mean) o la desviación estándar (sd). Con el comando data.frame() se puede crear una tabla de datos con los promedios por columna (poblaciones). Recuerde: si necesita la información de otro estrato, por ejemplo estadísticas a nivel de Sitio y no de Región, es necesario utilizar el comando setPop(x) <- region Antes de ejecutar el comando basic.stats y los demás comandos para obtener promedios.
Ho.prom <- apply(bs$Ho, MARGIN = 2, FUN = mean) #obtiene el promedio por población, de la heterocigosidad observada.
Ho.sd <- apply(bs$Ho, MARGIN = 2, FUN = sd) #obtiene la desv estandar de la heterocigosidad observada.
#Crear la tabla de datos con promedios se utiliza:
data.frame(Ho.prom, Ho.sd) Ho.prom Ho.sd
BRUNCA 0.498455 0.1404518
CHOROTEGA 0.540275 0.1271496Este resultado nos muestra que la heterocigosidad en Chorotega es levemente mayor (aproximadamente un 2% mayor). No obstante, la desviación estándar calculada entre loci, muestra que hay mucha diversidad entre loci, por lo que es difícil concluir que hay diferencias estadísticas entre las regiones. La conclusión más correcta es que la diversidad genética para el maíz criollo entre ambas regiones es similar. No hay pruebas estadísticas paramétricas para comparar estimativas de diversidad entre regiones, ya que no tienen una distribución estadística conocida. Usualmente se utilizan métodos de bootstrap para comparar, los cuales veremos más adelante.
Recuerde que al trabajar con objetos que son listas, se puede acceder a cada elemento de la lista con el símbolo de dólar $, de la manera objeto$elemento. En el ejemplo anterior, el objeto es bs y el elemento es Ho, por lo que accede solo a la heterocigocidad observada por población y por locus. También se puede obtener Fis, Hs, entre otros al usar el símbolo $, por ejemplo Hs.prom <- apply(bs$Hs, MARGIN = 2, FUN = mean) y fis.prom <- apply(bs$fis, MARGIN = 2, FUN = mean). Recuerde que cada vez que se presenta un promedio este debe venir acompañado con su desviación estándar o un intervalo de confianza. La desviación estándar puede ser calculada con el mismo comando, al cambiar el argumento FUN = sd.
3.4. Intervalos de confianza coeficiente endogamia
Frecuentemente, en un artículo científico se tiene como objetivo estimar la endogamia de la población o si se considera necesario para justificar ciertos resultados, es recomendable presentar el valor junto con un intervalo de confianza estadístico, con el fin de determinar si la estimación es significativamente distinta de 0 o para comparar si es distinta a la de otra población dada. Para ello, se puede usar la función boot.ppfis(x, nrep = 999). Esta función realiza un muestreo de 999 repeticiones por bootstrap a nivel de loci, con el fin de generar una distribución de valores de endogamia, a la cual se le extrae el cuantil 2.5% y el 97.5%, correspondientes al límite inferior y al superior de un intervalo de confianza al 95%.
Si el intervalo entre el límite inferior y el superior no incluye al 0, se considera significativamente distinto a 0. De igual forma, si esos intervalos no se traslapan entre las distintas poblaciones, se puede afirmar que son significativamente distintas. A continuación, se calcula el coeficiente de endogamia promedio para cada región como se explicó previamente y su intervalo de confianza 95%.
Fis.prom <- apply(bs$Fis, MARGIN = 2, FUN = mean) #obtiene el promedio por población, del coeficiente de endogamia.
#Crear la tabla de datos con promedios se utiliza:
data.frame(Fis.prom) Fis.prom
BRUNCA 0.375715
CHOROTEGA 0.332915boot.ppfis(magen,nboot = 999)$call
boot.ppfis(dat = magen, nboot = 999)
$fis.ci
ll hl
1 0.3072 0.4491
2 0.2824 0.3995A partir del resultado anterior, podemos determinar que la región Brunca posee un coefiente de endogamia significativo promedio de 0.376, con un intervalo 95% entre \left[0.303 ; 0.450\right]. Por otro lado, Chorotega también posee un coeficiente de endogamia significativo, con un coeficiente promedio de 0.333, con intervalo 95% entre \left[0.276; 0.402\right]. Sin embargo, no existen diferencias significativas entre los coeficientes de endogamia entre ambas regiones, ya que sus intervalos de confianza al 95% traslapan.
3.5. Número de alelos efectivos (Ae)
Esta estimación de diversidad genética, a diferencia del número de alelos, considera la frecuencia de cada alelo. Se estima para cada locus como \frac{1}{\sum{p_i}} donde p_i es la frecuencia de cada alelo. Dentro de las librerías de R no existe un comando que calcule Ae, por lo que debe realizarse el cálculo de forma manual. Para esto, se necesita la frecuencia de cada alelo por población, se obtienen con el comando pop.freq(x), donde x es el objeto genind. Con las frecuencias, se utiliza el comando lapply, que va a aplicar una función codificada por nosotros que llamaremos Ae. Esta utiliza la fórmula descrita para calcular el número de alelos efectivos.
alfrec <- pop.freq(magen) #crear objeto alfrec con frecuencias alélicas
Ae <- function (x){ 1/sum(x^2)} #crear la función con la fórmula para estimar Ae
#utilizar lapply en la lista de frecuencias y aplicar la función creada Ae
out <- lapply(alfrec, function(X) {Y <- as.data.frame(t(apply(X, MARGIN=2,FUN= Ae)))})
#Crear tabla de datos con marcadores en filas y columnas poblaciones
res.Ae <- do.call(rbind, out) # rbind permite unir filas.
res.Ae BRUNCA CHOROTEGA
bnlg1046 8.346026 6.572575
bnlg1191 12.482000 15.184430
bnlg1194 19.169065 21.400139
bnlg1265 14.231579 15.868545
bnlg1449 7.396253 7.580052
bnlg1523 9.515204 8.467183
bnlg1740 6.068834 11.026968
bnlg1890 6.981067 7.636089
bnlg1917 4.448959 5.448358
bnlg2305 8.142751 7.686078
phi015 3.775163 3.960958
phi031 3.197343 2.788605
phi064 9.224199 8.449608
phi072 3.376249 3.912630
phi078 3.244880 3.976607
phi089 2.045709 1.857790
phi093 4.678817 4.146252
phi109188 3.464335 4.847074
phi116 3.786604 5.246238
phi96100 3.416066 4.819662#para obtener el promedio por población (columnas), se utiliza:
colMeans(res.Ae) BRUNCA CHOROTEGA
6.849555 7.543792 Al igual que en resultados anteriores vemos que el número efectivo de alelos es levemente mayor en Chorotega. De nuevo, la diferencia parece ser leve. Un aspecto muy interesante es que el número de alelos promedio calculado para cada región fue superior a 17, pero el número efectivo de alelos es cercano a 7. Esto quiere decir que hay varios alelos (aproximadamente 7) con frecuencias altas y muchos más alelos de frecuencias bajas. En promedio, las regiones Brunca y Chorotega se comportan como regiones con 6 a 7 alelos con frecuencias similares. Siete alelos sigue siendo un número alto, en especial cuando se trabaja con más de 8 marcadores o loci. Tener menos de 4 alelos efectivos, usualmente se considera bajo, ya que no se tiene mucha información por locus.
3.6. Alelos privados
Son aquellos alelos que son exclusivos de una u otra población. Se utilizan más que todo para determinar diferencias entre poblaciones en estudio o buscar loci que puedan servir para clasificar individuos en una u otra población. Se obtienen con la función private_alleles (x, form = locus~. , count.alleles = F, report = "table"). Se debe indicar el objeto genind (x), si se quiere por locus (form = locus~.) o por alelo (form = alleles~. es la opción por default). Con count.alleles = F, cada alelo privado se contará una vez, independientemente de la dosis, es decir, en lugar del número observado de alelos privados. Con report = "table", muestra los resultados en una tabla.
#alelos privados por alelo
pal <- private_alleles(magen, count.alleles = F, report = "table")
pal bnlg1046.213 bnlg1046.215 bnlg1046.205 bnlg1046.217 bnlg1046.177
BRUNCA 1 0 0 0 0
CHOROTEGA 0 1 1 1 1
bnlg1046.208 bnlg1046.183 bnlg1046.211 bnlg1191.219 bnlg1191.231
BRUNCA 0 0 0 1 1
CHOROTEGA 1 1 1 0 0
bnlg1191.179 bnlg1194.215 bnlg1194.158 bnlg1194.219 bnlg1194.154
BRUNCA 0 1 1 1 1
CHOROTEGA 1 0 0 0 0
bnlg1194.223 bnlg1194.152 bnlg1194.201 bnlg1194.224 bnlg1194.166
BRUNCA 1 0 0 0 0
CHOROTEGA 0 1 1 1 1
bnlg1194.216 bnlg1194.182 bnlg1265.198 bnlg1265.239 bnlg1265.246
BRUNCA 0 0 1 1 1
CHOROTEGA 1 1 0 0 0
bnlg1265.212 bnlg1265.202 bnlg1265.208 bnlg1265.231 bnlg1449.163
BRUNCA 0 0 0 0 1
CHOROTEGA 1 1 1 1 0
bnlg1449.166 bnlg1449.135 bnlg1449.133 bnlg1449.138 bnlg1523.201
BRUNCA 0 0 0 0 1
CHOROTEGA 1 1 1 1 0
bnlg1523.139 bnlg1523.141 bnlg1523.265 bnlg1523.269 bnlg1523.235
BRUNCA 1 1 0 0 0
CHOROTEGA 0 0 1 1 1
bnlg1523.260 bnlg1523.233 bnlg1523.196 bnlg1523.231 bnlg1523.194
BRUNCA 0 0 0 0 0
CHOROTEGA 1 1 1 1 1
bnlg1523.250 bnlg1523.198 bnlg1740.150 bnlg1740.189 bnlg1740.168
BRUNCA 0 0 1 1 0
CHOROTEGA 1 1 0 0 1
bnlg1740.135 bnlg1740.166 bnlg1740.160 bnlg1740.130 bnlg1740.144
BRUNCA 0 0 0 0 0
CHOROTEGA 1 1 1 1 1
bnlg1740.179 bnlg1740.138 bnlg1890.204 bnlg1890.189 bnlg1890.197
BRUNCA 0 0 1 1 1
CHOROTEGA 1 1 0 0 0
bnlg1890.142 bnlg1890.132 bnlg1890.214 bnlg1890.140 bnlg1890.130
BRUNCA 1 0 0 0 0
CHOROTEGA 0 1 1 1 1
bnlg1890.163 bnlg1890.136 bnlg1890.165 bnlg1890.147 bnlg1890.145
BRUNCA 0 0 0 0 0
CHOROTEGA 1 1 1 1 1
bnlg1917.176 bnlg1917.164 bnlg1917.192 bnlg1917.136 bnlg1917.184
BRUNCA 1 1 1 0 0
CHOROTEGA 0 0 0 1 1
bnlg1917.88 bnlg1917.133 bnlg1917.145 bnlg1917.139 bnlg1917.166
BRUNCA 0 0 0 0 0
CHOROTEGA 1 1 1 1 1
bnlg1917.168 bnlg1917.150 bnlg1917.180 bnlg2305.206 bnlg2305.185
BRUNCA 0 0 0 1 1
CHOROTEGA 1 1 1 0 0
bnlg2305.200 bnlg2305.204 phi015.100 phi064.112 phi064.88 phi064.90
BRUNCA 0 0 0 0 0 0
CHOROTEGA 1 1 1 1 1 1
phi064.108 phi072.166 phi072.146 phi078.126 phi089.89 phi093.299
BRUNCA 0 0 0 0 1 0
CHOROTEGA 1 1 1 1 0 1
phi093.294 phi109188.151 phi109188.161 phi109188.174 phi116.177
BRUNCA 0 1 0 0 0
CHOROTEGA 1 0 1 1 1
phi116.169 phi96100.285 phi96100.287
BRUNCA 0 0 0
CHOROTEGA 1 1 1Para caluclar el total de alelos privados, debemos sumar por filas, y esto lo podemos hacer con el commando apply().
apply(pal,1,sum) #promedio de alelos privados por población BRUNCA CHOROTEGA
28 77 A partir de esto, se puede concluir que la región Brunca tiene 28 alelos exclusivos de esa población, mientras que en la región Chorotega, existen 77 alelos exclusivos de esa región.
Finalmente podemos calcular el numero de alelos privados por locus,
##alelos privados por locus (por SSR)
paloc<-private_alleles(magen, form = locus~. , count.alleles = F, report = "table")
paloc bnlg1046 bnlg1191 bnlg1194 bnlg1265 bnlg1449 bnlg1523 bnlg1740
BRUNCA 1 2 5 3 1 3 2
CHOROTEGA 7 1 6 4 4 10 8
bnlg1890 bnlg1917 bnlg2305 phi015 phi064 phi072 phi078 phi089 phi093
BRUNCA 4 3 2 0 0 0 0 1 0
CHOROTEGA 9 10 2 1 4 2 1 0 2
phi109188 phi116 phi96100
BRUNCA 1 0 0
CHOROTEGA 2 2 2Aquí podemos ver que el locus bnlg1194 tiene 5 alelos privados.
Referencias consultadas
Bedoya, C. (2012). Estudios de diversidad genética en poblaciones de maíz (Zea mays L.) evaluadas con microsatélites. Tesis doctoral. Universidad de las Islas Baleares, España. 223p.
Fuchs, E. & Barrantes, G. (2015). El lenguaje estadístico R aplicado a las ciencias biológicas. Editorial UCR, San José, Costa Rica. 197p.
Jiménez, R. (2014). Caracterización de las razas criollas e indígenas de maíz colombiano por medio de marcadores moleculares SSR. Tesis para optar por el grado de Magister en Ciencias Biológicas. Universidad Nacional de Colombia Facultad de Ciencias Agropecuarias. 70p.
Grünwald, N., Kamvar, Z. & Everhart, S. (2015). Population genetics in R. Corvallis, Oregon, USA. Recuperado de http://grunwaldlab.github.io/Population_Genetics_in_R/Data_Preparation.html [Consulta 02 feberero 2017]
Paradis E. 2010. pegas: an R package for population genetics with an integrated-modular approach. Bioinformatics 26: 419-420.
Peakall, R. & Smouse, P.E. (2012). GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research-an update. Bioinformatics 28(1): 2537–2539. R Core Team (2013).
R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.